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Abstract 

The ability of the three-dimensional Navier-Stokes 
method, PAB3D, to simulate the effect of Reynolds num- 
ber variation using non-linear explicit algebraic Reynolds 
stress turbulence modeling was assessed. Subsonic flat plate 
boundary-layer flow parameters such as normalized veloc- 
ity distributions, local and average skin friction, and shape 
factor were compared with DNS calculations and classical 
theory at various local Reynolds numbers up to 180 million. 
Additionally, surface pressure coefficient distributions and 
integrated drag predictions on an axisyinmetric nozzle af- 
terbody were compared with experimental data from 10 to 
130 million Reynolds number. The high Reynolds data was 
obtained from the NASA Langley 0.3m Transonic Cryogenic 
Tunnel. There was generally good agreement of surface 
static pressure coefficients between the CFD and measure- 
ment. The change in pressure coefficient distributions with 
varying Reynolds number was similar to the experimental 
data trends, though slightly over-predicting the effect. The 
computational sensitivity of viscous modeling and turbu- 
lence modeling are shown. Integrated afterbody pressure 
drag was typically slightly lower than the experimental data. 
The change in afterbody pressure drag with Reynolds num- 
ber was small both experimentally and computationally, 
even though the shape of the distribution was somewhat 
modified with Reynolds number. 

Introduction 

Current focused program efforts are considering Reynolds 
number scaling a significant aspect of aircraft testing and 
development. Wing aerodynamics and flow about propul- 
sion systems can have considerable sensitivity to varying 
Reynolds number. Most of the sub-scale wind tunnel testing 
occurs at Reynolds numbers below that of flight conditions; 
therefore, the ability of computational fluid dynamics (CFD) 
to simulate higher Reynolds number flow is of importance. 

Previous to the development of cryogenic test techniques 
for achieving high Reynolds numbers in wind tunnel facili- 
ties, little fundamental research data had been available for 
the evaluation of any theoretical methods to predict these 
effects. Several years ago. during the development phase of 
cryogenic testing techniques at the NASA Langley Research 
Center; two sets of simple axisy mine trie nacelle models were 
built and tested in what, was then known as the l/3m Pilot 
Transonic Cryogenic Tunnel (now the 0.3m Transonic Cryo- 
genic Tunnel). This was some of the first set. of test data 


* Senior Scientist. Component Integration Branch. Aerodynamics Division. Senior 
Member AIAA. 

Copyright (?) 199G by the American Institute of Aeronautics and Astronautics. Inc. 
No copyright is asserted in the United States under Title 17. U S. Code. The U.S. 
Government has royalty-free license to exercise all rights under the copyright claimed 
herein for Governmental purposes. All other rights are reserved by the copyright owner. 


for nozzle-boat tail geometries taken over a large range of 
Reynolds numbers, refs. 14. 

The current investigation assesses the capability of the 
Navier-Stokes method PAB3D, version I3S, (refs. 5 8) using 
non-linear algebraic Reynolds stress turbulence models to 
predict the Reynolds number effects on the flow about a 
nozzle boattail. and simulate a 5 meter flat plate at very 
high Reynolds numbers. Comparisons were made with wind 
tunnel data for the boattail geometry and boundary layer 
profiles, shape factor, and skin friction with DNS data and 
textbook equations for incompressible flat plate flow. 

Nomenclature 

A umx maximum body cross-sectional area. 

0.78539 in 2 

p 

Cpf pressure drag coefficient. ^ ~ q maj 

Cp average skin friction coefficient. 

bb £ r,,Ai 

Cj local skin friction coefficient. r tr /q^c 

Cp pressure coefficient, ^ 

d m body maximum diameter. 1.0 in. 

F axial force along body axis 

fp near-wall damping function for linear 

K-e 

GS Gatski-Speziale 

H \ 2 boundary layer shape factor. 61/62 

Hyi boundary layer shape factor. 63/62 

h\ physical height of first computational 

grid from a wall 

K turbulent kinetic energy 

l integration length of flat plate 

L model reference length 

M Mach number 

Reynolds number based on model refer- 
ence length 

n direction normal to wall 

P production term for turbulent kinetic 

energy 

p static pressure. Pa 

q dynamic pressure. Pa 

Rl Reynolds number based on flat plate 

integration length 

Rj cell turbulent Reynolds number. K 2 /vt 

R v Reynolds number based on distance; x. 

u^exjv 
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R^ displacement thickness Reynolds number, 

UocSl/l' 

Rq momentum thickness Reynolds number, 

t time 

S strain tensor 

SZL Shih, Zhu & Lumley 

U magnitude of local velocity, ti* 2 

u stream-wise velocity 

Ufc cartesian velocity components 

u + law-of- the- wall coordinate, uju T 

friction velocity, \Z(r w /p) 
t/t/"*" nondimensional shear stress, uV/ti r ^ 

W vorticity tensor 

x stream-wise distance 

law-of-the-wall coordinate. nu T /v 
z vertical distance 

Al incremental distance on flat plate 

61 boundary layer displacement thickness 

62 boundary layer momentum thickness 

63 boundary layer energy thickness 

e turbulent dissipation 

p laminar viscosity 

fit turbulent viscosity 

fi w local laminar viscosity at the wall 

v kinematic viscosity, fij p 

p density 

r shear stress 


Computational Procedure 
Governing Equations 

The code used was the general three dimensional (3-D) 
Navier-Stokes method PAB3D, version 13S. This code has 
several computational schemes, different turbulence models, 
and viscous stress models that can be utilized, as described 
in more detail in refs. 5 through 8. The governing equations 
are the Reynolds- averaged simplified Navier-Stokes equa- 
tions (RANS) obtained by neglecting all stream-wise deriva- 
tives of the viscous terms. The resulting equations are 
written in generalized coordinates and conservative form. 
Viscous model options include k-thin layer, j-thin layer, 
jk-uncoupled and jk-coupled simulations. Typically the thin- 
layer viscous assumption of the full 3-D viscous stresses is 
utilized. Experiments such as the investigation of super- 
sonic flow in a square duct was found to require fully cou- 
pled 2 directional viscosity to properly resolve the physics 
of the secondary cross-flow. The Roe upwind scheme with 
first, second, or third order accuracy can be used in evalu- 
ating the explicit part of the governing equations and the 
van Leer scheme is used to construct the implicit operator. 
The diffusion terms are centrally differenced and the inviscid 
flux terms are upwind differenced. Two finite volume flux- 
splitting schemes are used to construct the convective flux 
terms. 

All solutions were developed using third-order accurate 
schemes for the convective terms, and second-order for the 
viscous diffusion terms, denoted by the first 3 in the nomen- 
clature in the figures and tables and the min-mod solution 
limiter, denoted by the second 2 in the nomenclature. Only 
the viscous model is varied in this study, denoted by the 
third number in the nomenclature. For completeness, a ta- 
ble of nomenclature designating the order of scheme, limiter, 
and viscous modeling is given below. Other solution limiters 
include van Albeda, Spekreijse-Venkat (S-V) and a modified 
S-V (ref. 9). Solution limiters influence solution convergence 
and final results. In some instances, such as a jet-plume sim- 
ulation, the van Albeda solution limiter is required to obtain 
a smooth converged solution. 


4> angular location of pressure orifices, deg 

Superscripts 
L laminar 

T turbulent 


Subcripts 

0 

CL 

fP 

l 

n 

t 

to 

sf 

w, wall 
00 


nozzle boattail component contribu- 
tion 

centerline 
flat plate 
laminar 

non-linear component 
turbulent 

free stream total condition 
skin friction contribution 
condition at the wall surface 
free stream condition 


Nomenclature 

Solution Limiter 

Viscous model 

311 

van Albeda 

k-thin layer 

312 

van Albeda 

jk-coupled 

313 

van Albeda 

jk-uncoupled 

321 

min-mod 

k-thin layer 

322 

min-mod 

jk-coupled 

323 

min-mod 

jk-uncoupled 

331 

S-V 

k-thin layer 

332 

S-V 

jk-coupled 

333 

S-V 

jk-uncoupled 

341 

modified S-V 

k-thin layer 

342 

modified S-V 

jk-coupled 

343 

modified S-V 

jk-uncoupled 


The code can utilize either a 2-factor or 3- factor numerical 
scheme to solve the flow equations. The 2- factor scheme is 
typically used as it requires 10 to 15 percent less memory 
as compared to the 3-factor scheme. The memory difference 
is dependent on the size of cross- planes of the specific grid 
being used. When the 2-factor scheme is used the orientation 
of the grid and predominate flow direction typically along 
the i grid index, such that the Roe scheme is utilized to 
sweep stream-wise through the computational domain and 
the van Leer scheme for the solution of the cross- plane 
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(i.e.. i = constant) of a 3-D problem. However solving a 
single-cell wide two-dimensional (2-D) mesh defined with 
the i direction of the grid oriented in the conventional 
stream-wise direction will typically converge slower using 
the Roe relaxation solution scheme compared to solving the 
equivalent problem with the van Leer scheme. Therefore the 
i and j directions of a 2-D mesh are swapped allowing the 
entire flow- field to be solved implicitly with each iteration. 
The explicit sweep is not used since only one cell exists in the 
i direction. The implicit scheme usually has a much higher 
rate of convergence and typically provides a solution using 
less computational time. 

Turbulence Simulation 

The turbulence model equations are uncoupled from the 
RANS equations and are solved with a different time step, 
typically 1/2. than that of the principle flow solution. A 
considerably lower principle Courant-Friedrichs-Levy (CFL) 
number is typically required to solve problems if both the 
main flow equations and turbulence equations are solved it- 
eratively using identical time rates. Larger time step differ- 
ences, e.g., 1/4 to 1/8, slow solution convergence further but 
result in identical final solutions. Flow solution transients at 
times require the turbulence equations time step to be re- 
ducted temporarily. Turbulence simulations are resolved at 
all grid levels, not just at the finest grid level. 

Version 13S of the PAB3D code used in this study has 
options for several algebraic Reynolds stress (ASM) turbu- 
lence simulations. The Standard model coefficients of the 
K — e equations were used as the basis for all the linear and 
non-linear turbulent simulations, ref. 10. Additionally, it is 
known that the eddy viscosity models produce inaccurate 
normal Reynolds stresses. Flat plate flow, as well as other 
more complex aerodynamic flows, are anisotropic. 

Successful implementation of the algebraic Reynolds 
stress models required the solution methodology for turbu- 
lent production term P of the underlying linear turbulence 
calculations to be modified. P depends on high order deriva- 
tives of the turbulent Reynolds stresses. Proper represen- 
tation of the stresses should be provided by face centered 
values, rather than the cell centered values. Previous at- 
tempts to implement non-linear turbulence models in the 
context of a cell centered eddy viscosity model worked only 
for 2-D problems and was unable to resolve 3-D flows. 

Linear K - £ equations The transport equation for the 
turbulent kinetic-energy. K. and the dissipation rate are 
written as: 


where P = and (C £ ~i=1.44, CL 2=1.92, C^^O.090). 

The damping function of Launder & Sharma. ref. 11, 
f fi = exp(^ -3.4 1/(1 -f Rt! 50. ) 2 ^. determined the behavior 
of e near the wall as a function of turbulent Reynolds number 
Rj = K} jus. The boundary conditions for e and K at 

the wall are £ wa n = and K tV(l u = 0. The 

stress components in linear turbulence models are developed 
with laminar and turbulent components. r t j = rjj 4- tJj. A 
generalization of Boussinesq's hypothesis redefines laminar 
and turbulent components are as follows: 

4 = - 2 ii L S u (3) 


where 


A L = -n L s kk 


and 



( du < . du A 

\dxj dxj ) 


(4) 


T 

The turbulent component, of the stresses r f j is repre- 
sented by the sum of linear (T/) and non-linear (T n ) compo- 
nents. The linear stress is rjj = A T 6jj — 2 p T S tJ where 

A t = \(pK -I- p T S kk ). The non-linear component of the 
turbulent, stresses are addressed in the following section. 

Non-Linear Turbulent Stress Equations- Three theories 
of explicit algebraic Reynolds stress models were imple- 
mented. The Reynold's stress contribution rjp used by 
Shih. Zhu, & Lumley (SZL), (ref. 12) is: 

T 5" = 20^- ( W ik S kj - S ik W kj ) (5) 

Gatski & Spcziale (GS), (ref. 13); 

rff = c;~ [0i(w, k s kj - s lk w kj ) + Ih(s ik s kj 

— (®) 
and Girirnaji (G). (ref. 14); 


T„ _ 




= 2C,:^r [-<? 2 (W, k S kj - S ik W kj ) + G 3 (S ik S kj 


de 

di 


de d (, L _ K 2 . de \ 

■**w k -or k { lv 


, r iZ-r — 

+ Crl K 12 K 




( 1 ) 


— ~ Sm n Sm n &i j ) 


where 


W; 


_ 1 / du } duj 
2 \dxj dx t 


(7) 


The convective terms are solved using third-order differ- 
encing. The diffusion terms are solved using second-order 
central differencing. 


S}j — j ~S kk 6jj 


The turbulent viscosity. /i . is defined as 


dK dK 
dt + dx k 


= P~£ + 


dx k 


£ dx k 


( 2 ) 
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where C £ = for solutions solving linear turbulence sim- 

ulations and equal to the variable function C J = /( S, W 1 K, e) 
for solutions involving algebraic Reynolds stress simulations. 
Functions for C J take the following forms for each of the 

ASM. 

Shih, Zhu & Lumley, (ref. 12): 

c . = i /( 6 . 5 + A ;^) (9 ) 

Gatski-Speziale, (ref. 13): 

C* = const. * {1 + X 2 )/{3 + X 2 + &x 2 i> 2 + c xp 2 ) (10) 

Ali U* . x« and fp are all different functions of the strain 
and vorticity tensors and are detailed in the references. 

Girimaji, (ref. 14): 


Solution Process 

Turbulent flow solutions using ASM and two-equation 
linear K — e model requires 23 words per grid point. The 
code speed is dependent on the turbulence model, thin-layer 
assumptions and numerical schemes. The following table 
are some options available in the code with C-90 timing in 
/iseconds/iteration/grid point. 


Solver Scheme 

Viscous 

Modeling 

Turbulence 

Modeling 

(3 r ^-order) 

Stress 

Center 

Timing C-90 
/is/ iter /grid 

2-factor 

j-k uncoupled 

Girimaji ASM 

Face 

23 

2- fact or 

k thin-layer 

Girimaji ASM 

Face 

20 

Diagonalization 

j-k uncoupled 

Girimaji ASM 

Face 

16 

Diagonalization 

k thin-layer 

Girimaji ASM 

Face 

14 

2-factor 

k thin-layer 

Gat ski k 
Speziale ASM 

Face 

19 

2-factor 

k thin-layer 

SZL ASM 

Face 

20 

2- factor 

k thin-layer 

Linear- Isotropic 

Face 

18 

2-factor 

k thin-layer 

Linear- Isotropic 

Cell 

17 

Diagonalization 

k thin-layer 

Linear- Isotropic 

Face 

12 


Gi - 


L\L 2 / [|£») 2 + 2rfe(I„) 2 ] 

L\UJ [(£®) 2 + |r/i(L 3 ) 2 + 2 % (£ 4 ) 2 ] 

-f + c ° 8 i j) 

_ § + cos ( J + t) 


for 7/1 = 0 : 
for L\ — 0; 

for D > 0: 

for D < 0 and b < 0; 
for D < 0 and b > 0. 


(ii) 


The variable G \ utilized by Girimaji is equal to — CJ. 
A compilation of the parameters used in Girimaji *s model 
can be found in the Appendix. Additional information is in 
reference 14. 

The solution processes for wall-bounded flows were 
equally robust for each of the models. Previous results, not 
published here, show Gatski-Speziale requiring lower CFL 
numbers for the solution of free-shear flows. Obtaining con- 
verged solutions using Gatski's C* were found to be problem 
dependent. Girimaji s Gi function appears to be extremely 
well behaved permitting for fairly high CFL numbers to 
used. 

Turbulent Trip Equation*- The technique used for initial- 
izing the viscous flow transition from laminar to turbulent is 
placing K and e profiles at user-specified lines or planes in the 
flowfield. The line or plane of the specified trip area is sur- 
veyed for the maximum and minimum velocity and vorticity 
along that line and a shape function from 0 to 1 is created of 
the form F — (J — f nun )/{fmar ~ fmin) where f is a prod- 
net of the velocity and vorticity / = u|W|, \ W\ = 2 W '^ . 

The turbulent kinetic energy profile is then K = a U F . 
where q is a free parameter determining the magnitude of 
the impulse as a percent of local total velocity. U . The typi- 
cal value specified by the user, and used for this paper, is 2% 
(or a = 0.02). The e profile is developed from the assump- 
tion that production P is equal to the dissipation e equaling 
vSfjTjj*-. The result of the initialization is seen as a 
spike in the K field of the solution. This initial turbulent, 
profile develops as permitted by the local flow conditions. 


Several parameters were used to gauge solution conver- 
gence. Local skin friction, shape factor and solution residual 
were monitored for convergence of the flat plate solutions. 
Total afterbody drag, nozzle pressure drag, and solution 
residual were used to determine the solution status at the 
coarse (144), medium (122), and fine (111) grid levels of the 
axisymmetric afterbody. The 144 abbreviation means divide 
number of i-cells by 1, number of j-cells by 4 and the number 
of k-cells by 4. Afterbody drag variance of less than 0.50 per- 
cent for several hundred iterations was achieved for all test 
cases. 

The conservative patch interface package of Pao and 
Abdol-Hamid (ref. 7) enables the code to properly transmit 
information between mis-matched block interfaces. Integer- 
to-one interfaces are considered a subset of the arbitrary 
block interface and do not need to be specified as such to 
the patching code. The patching program is a preprocessor 
that writes a connectivity data base prior to the start of the 
first solution. Each entry to the patch data base contains 
cell face areas and indices relating that cell with all other 
cells that will share momentum flux information. The data 
base information is automatically re-allocated internal to the 
code during mesh sequencing. As a result, each block can 
be sequenced at different levels and the correct interface 
information is maintained at the cell level. However, it is 
important to note that features in the flow developed on one 
side of an interface should not be obliterated on the other 
side due to an excessive grid density mis-match. 

Third-order continuity in transmitting the fluxes across 
block boundaries is maintained by the code: lower order 
continuity may be specified by the user if required. As 
with most Navier-Stokes methods of the typo, equal cell size 
spacing on either side of an interface in directions normal 
to the interface should be maintained regardless of the mesh 
sequencing level of the block. 

Boundary Conditions 

For this study, solid walls were treated as no-slip adiabatic 
surfaces. The solid wall boundary condition was satisfied 
by setting the momentum flux of the solid wall cell face 
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to zero. A boundary condition for the Riemann invariants 
along the characteristics was specified for the free-st.ream 
inflow face and the lateral free-stream outer boundary of 
the flow domain. An extrapolation boundary condition was 
applied on the downstream outflow face. The axisymrnetric 
flow assumption for the single-cell grids was implemented 
by placing flow symmetry conditions to the lateral side 
boundaries of the computational domain. 

Results and Discussion 
Subsonic Flat Plate 

Flat Plate Grid The 5 m flat, plate niultiblock grid had 
an H-type mesh topology, with the blocking sketched in 
figure 1. The computational domain included inflow block 
extending 1 meter upstream from the leading edge of the 5 
m flat plate. The initial stream- wise grid spacing at the 
leading edge of the plate was 1.x 10 4 m and was exponentially 
stretched from the leading edge to the trailing edge at a 
rate of 5% with a total of 1C1 grid points. The first cell 
height was 1.0 x 10^ m fixed at both ends of the plate 
and exponentially stretched from the surface to the outer 
boundary at a rate of 11% with a total of 121 grid points. 
The upper boundary was 2 m away and the lateral width 
of the grid of 0.01 rn. All three blocks had dimensions 
of 81 x 121. Tripping to turbulent flow simulation occurred 
around R x = -3 million or R^\ = 900. corresponding to a 
physical distance of approximately 9 mm downstream of 
the plate leading edge. This allowed for laminar flow to 
occur over roughly 32 computational cells before tripping to 
turbulent flow. Grid cell counts were divisible by four to 
allow a minimum of 2 levels of grid sequencing. 

Boundary Layer Characteristics — Figure 2 shows the 
Reynolds number based on length variation with distance 
from the leading edge. The Reynolds number at the plate 
trailing edge was approximately 180 million. Note that the 
plot is a log-log type with the symbols indicating the stream- 
wise distribution of the grid points. The high Reynolds 
number was obtained through increasing the free-stream 
total pressure, rather than physically lengthening the fiat 
plate geometry. The normalized velocity and shear stress 
distributions at R# = 1420 and 100.000 are shown in fig- 
ures 3 and 4. The comparisons at R$ = 1420 are com- 
pared with the DNS calculations of Spalart.. ref. 15. and 
at R# = 100,000 are compared with the classical flat plate 
equations. All three ASM match fairly closely the DNS 
calculation shown in figure 3, with the Girimaji model fol- 
lowing the closest, in the buffer region. All three models 
were slightly above the DNS at. the edge of the bound- 
ary layer. Similarly. Girimaji best fit the DNS stress pro- 
file. u'V + = (du/dzjCfJuK^/e/ur. though all three ASM 
were generally a good match. The high Reynolds num- 
ber comparisons, figure 4. at R 9 = 100,000. approximately 
TV/?, = 90 million, have trends fairly consistent with the clas- 
sical flat plate boundary layer flow equations. The stress 
profiles, figure 4(1)). have similar lower level behavior (be- 
low = 50) as the lower Reynolds number profiles and a 
greatly flattened region of constant, stress below the bound- 
ary layer edge around — 30,000. The grid had typically 2 
cells less than = 2.5 and about 3G cells in the boundary 
layer at Ry = 1420. 

Flat Plate Skin Friction Figures 5 and G are a compari- 
son of classical flat plate theories for local and average skin 


friction with the three ASM solutions. The equations for the 
local skin friction comparisons were: 


c f 


( 0.G64A/ET, 

{ 0.0590 R x ~h, 


{ 0.455/in 2 (O.OGRx) > 
The equations for the average 


B Iasi us ; 

^ th power law; 

White- “Exact" theory, 
skin friction were: 


( 12 ) 


Cf = 


’ 1.328/ v/R^, 

QAbbl(log\™(R L )-A!R L ), 

* _ l 

0.074 R l S - A/R l , 

. 0.523//n 2 (0.00fli). 


Blasius: 

Transition; 

ith power law; 

White- “Exact" theory. 

(13) 


where A = Rentier, ~^F { )- Cf) = 1-328/ y/Rcrit • Cf t = 
0.074 [Rcrit)-*. 

R cr it i s the local Reynolds number at the point of transi- 
tion from laminar to turbulent flow. Transition was defined 
as the point at which the shape factor H\ 2 first fell below 
2.3. Local skin friction and average skin friction coefficients 
and normalized turbulent viscosity are plotted in figures 7, 8 
and 9. respectively for all three of the algebraic Reynolds 
stress models. Girimaji. SZL and GS ASMx predict sim- 
ilar and consistent skin friction characteristics throughout, 
the Reynolds number range. All three models were virtually 
identical in local and average skin friction for the laminar 
flow that developed upstream of the transition trip point at 
R x = 300,000. Downstream of the trip, the Girimaji model 
developed slight higher local skin friction that the other two 
ASM. with subsequently higher average skin friction. All 
three models departed from the 1 /5th power theory for local 
skin friction at Reynolds numbers above 20 million. The skin 
friction predicted by Girimaji's model was slight ly above the 
higher Reynolds number theory of White, while the other 
two tracked slightly low. 


The trend of average skin friction through transition 
to turbulent flow was similar between the three models 
and followed the 1 /5th power theory very closely until, 
again departing around 20 to 30 million Reynolds number, 
figure G. Figure 7 is a plot of the growth of turbulent 
viscosity normalized by the local laminar viscosity with R 
Girimaji’s model predicts the highest level of normalized 
turbulent viscosity, though all three models are very similar 
in level and rate of growth. 


Boundary Layer Shape Factors - All three ASM have very 
similar shape factor H \ 2 trends as shown in figure 8. The 
first. 8 or so computational cells were neither laminar nor 
turbulent as the solution developed. The subsequent 28 
cells matched the theoretical laminar characteristics very 
closely. The theoretical turbulent shape factor was not 
closely achieved until around R x — 20 million. Even though 
transition from laminar flow occured relatively quickly, for- 
mation of a turbulent shape factor dost’ to the theoretical 
shape required some distance to achieve. All three models 
very closely match the turbulent shape factor of H\o = 1-27 
at very high Reynolds numbers. 


Overall, all three non-linear turbulence models appear to 
be consistent and well behaved turbulent, flat, plate proper- 
ties up to Reynolds numbers of 180 million. 
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Axisymmetric Afterbody 

Test Facility— The second test case was an axisymmet- 
ric geometry that was part of a series of models tested 
in both the Langley l/3m Pilot Transonic Cryogenic Tun- 
nel and the 16-Foot Transonic Tunnel. The Pilot Tunnel 
had an octagonal test section with slots at the corners of 
the octagon and is essentially a scale model of the Lang- 
ley 16-Foot Transonic Tunnel test section, ref. 16. The test 
medium for the cryogenic tunnel was nitrogen cooled by liq- 
uid nitrogen. High Reynolds number data were obtained in 
the 0 . 3 m tunnel through a combination of cryogenic free- 
stream temperatures and free-stream total pressure that are 
independently controllable. Approximately 5 atm. of pres- 
sure and 100 K total temperature produced a unit Reynolds 
number of 260 million/meter. 

The experiment was conducted over a range of tempera- 
tures from approximately 100 K to 300K and pressures from 1 
to 5 times the standard atmospheric level. Several settings 
of free-stream total temperatures or pressures can result in 
identical settings of Reynolds number. Surface pressure co- 
efficients and nozzle boattail drag were shown to be simi- 
lar regardless of the temperature/pressure combinations that 
created equivalent Reynolds numbers, ref. 2. High Reynolds 
number simulations with the CFD method were obtained 
through increased total pressure rather than through a com- 
bination of free-stream total pressure and cryogenic tem- 
peratures. Though data were obtained over range of Mach 
number from 0.6 to 0.9, only the M = 0.9 data is compared 
with the CFD in this paper. The following is a table of 
conditions for experimental data obtained at M = 0.9 for 
the L/d m = 16.0 model. One atmosphere is defined at 
0.101325 MPa (14.703 psi). 


A/oc 

r (0 ,K(R) 

p* 0 .atm 

■V/!. X 10 6 

.903 

106 (191) 

4.98 

128 

.908 

118 (212) 

3.98 

87 

.901 

119 (214) 

2.9$ 

64 

.911 

118 (212) 

2.47 

55 

.910 

118 (212) 

1.97 

43 

.904 

119 (214) 

1.49 

32 

.903 

118 (212) 

1.24 

27 

.899 

312 (562) 

4.97 

28 

.899 

308 (554) 

3.79 

22 

.902 

308 (554) 

2.48 

14 

.901 

307 (553) 

1.23 

7 


Geometry - The configuration used for this study was 
one of six models that were built for the original Reynolds 
number study, ref. 1 . Four models with differing boattail 
geometry were associated with a body length of 8 inches 
from the nose to the start, of the boattail (characteristic 
length) and two models with a characteristic length of 16 
inches. The boattail geometries had circular arc. circular 
arc-conic, or contoured profiles. This investigation utilized 
the circular arc with a leugth-to- maximum-diameter ratio 
(fineness ratio) of 0.8 boattail. Figure 9 is a photograph 
of the model mounted in the pilot tunnel. The nose of 
the model was a 28° cone 1.7956 inches long fairing to the 
cylindrical body via a 1.3615 inch radius circular arc whose 
center is 2.125 downstream of the model nose and 0.8615 
inches below the model centerline. The circular arc failing 
is tangent at its endpoints to the conical nose (1.7956 inches 
from the nose) and cylindrical body (2.125 inches from the 


nose). The mode! was sting mounted with the diameter 
of the sting being equal to the model base diameter. The 
length of the constant diameter portion of the sting (6.70 
inches measured from the nozzle connect station) was such 
that, based on the work of Cahn. ref. 17. there should be 
no effect of the sting flare downstream of the nozzle trailing 
edge on the boattail pressure distributions. 

The axisymmetric afterbody grid utilized H-0 type mesh 
topology with all block dimensions that were divisible by 4. 
The mesh was gridded with a single cell 5 degree wide wedge 
grid with the stream- wise flow direction oriented along the j 
index to utilize the implicit flow solver in the code for faster 
solution convergence. The body was described using 100 
cells extending from the leading edge of the nose to the nozzle 
connect station. There were 80 cells extending from the 
nozzle connect station to the nozzle boattail trailing edge. 

The free-stream conditions for axisymmetric CFD cases 
were M = 0.9, 7 / 0 = 540R using air at 7 = 1.4. The first 
cell height of each configuration's grid was different for 
each free-stream Reynolds number according to the following 
schedule. 


Reynolds number 

/^O. atm. (psi.) 

h \ (inches) 

7 

1.2 (17.8) 

6xl0~ 6 

55.2 .... 

9.52 (140.) 

SxlO -6 

128.3 .... 

22.1 (325.) 

2xl0 -6 


The wind tunnel models were constructed of cast alu- 
minum with stainless-steel pressure tubes cast as an integral 
part of the model. The model was instrumented with 30 
pressure orifices in three rows of 10 orifices each. The 1 inch 
diameter of the model physically precluded the placement 
of all 30 orifices along the same row. The following is a 
tabulation of the non-dimensional orifice locations. 


x/dm for L /d m = 16 at 

<£ = 0° 

d> = 120° 

<t> = 240° 

-0.4491 

-0.4660 

-0.4561 

-.1637 

-.2201 

-.1552 

-.0600 

-.1281 

-.0590 

.0337 

-.0260 

.0390 

.1268 

.0744 

.1342 

.2279 

.1729 

.2713 

.3210 

.26% 

.3718 

.4199 

.3679 

.4680 

.5231 

.4640 

.5749 

.6279 

.6758 

.7304 


Grui convergence Figures 10 and 11 show grid sensitiv- 
ity of the Girimaji ASM at M = 0.9 at the lowest and highest 
Reynolds number for this test case. Nr ( , =7 and 128 million, 
respectively. These sensitivities were relatively consistent for 
the other turbulence models and other viscous models inves- 
tigated. A few exceptions occurred where the coarse grid 
solution did not converge, but the following medium and 
fine grid solutions converged and the results were similar in 
nature as those shown in figures 10 and 11 . All solutions 
were fairly well grid converged and solution converged. Ini- 
tial inspection of figure 11 . the coarse grid solution has the 
closest match with the data. Further refinement of the grid 
revealed this solution to not be grid converged. Converged 
solutions for this geometry appear to require between 40 
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to 80 cells along the nozzle boattail to adequately predict 
the shock-separated flow reasonably accurately. 

Low Reynolds number Computations - Figures 12 through 
16 are low Reynolds number calculations showing the ef- 
fect of turbulence model, turbulent trip location and viscous 
model on pressure coefficient and turbulent kinetic energy 
distributions. 

In figure 12 all the calculations were performed with us- 
ing a single thin-layer viscous model, i.e.. k-thin layer for 
this mesh, the min-mod solution limiter, and a turbulent 
trip point, trip 1. approximately 0.031 inches (0.08 cm) 
downstream of the nose. The three ASM predicted a shock 
strength slightly weaker than the data and a pressure re- 
covery slightly lower than the data. The Standard K — e 
model, in this instance, appears to have better agreement 
with the data closely matching peak negative pressure and 
recovered to a static pressure only slightly above that of the 
data at the boattail trailing edge. Figure 13 is a plot of 
the peak turbulent kinetic energy for each turbulence model 
using the same parameters as the calculations in figure 12. 
For clarity, two areas of the axisymmetric body are de- 
tailed, the region downstream of the nose where the tur- 
bulent trip occurs and the region around nozzle boattail. 
The large spike in K/a 2 just downstream of x/d m = -16. 
is the turbulent trip impulse in k. None of the four tur- 
bulent models tested developed turbulent flow immediately 
downstream of the trip. The Standard K — e linear model 
developed turbulence first as seen by the rise in K/a 2 around 
x/d m = —15.7. The Girimaji and SZL ASM became turbu- 
lent, around x/d m = -15.3. and GS became turbulent the 
furthest downstream at xjd m — —14.4. 

Early studies simulating the incompressible flat plate 
flow displayed similar characteristics. If the turbulent 
trip was placed upstream of the critical flow point, tur- 
bulence would not develop immediately downstream of the 
trip. Conversely, turbulence would develop immediately 
when the turbulent trip was placed downstream of the 
critical flow point. Considering this, a different turbu- 
lent trip point was chosen roughly between the furthest 
upstream and downstream turbulent development, points 
noted previously and solutions re-developed for the three 
ASM. Figures 14(a) through (c) show that downstream 
of the cone-cylinder transition of body shape, approx. 
x/dm — — 13, despite the different initial development of 
turbulence, (trip 1 upstream trip® x/d m — —15.969 ver- 
sus trip 2 downstream tripC* x/d m = —15.000) no signif- 
icant changes occur in the peak turbulent kinetic energy. 
Figure 15 is representative of the lack of influence on 
static pressure coefficient distribution on the nozzle boat- 
tail between the two turbulent trip [joints using the min- 
mod solution limiter. Further parametric studies are needed 
to determine the boundary layer behavior using other solu- 
tion limiters with changes in the laminar-to-turbulont flow 
regions. 

Figure 16 is a study of the effect of different viscous 
models on the flow oil the afterbody. Three calculations 
were performed using k-thin layer (321): j-k viscosity coupled 
(322): and j-k viscosity uncoupled (323) viscosity models 
with Girimaji ASM at 7 million Reynolds number. The 
use of j-k viscosity appeal's to improve' the comparison with 
experimental data by creating a shock slightly stronger and 
further downstream than the k-thin layer calculation, in 


addition to slightly raising the pressure recovery in the 
region of separated flow. As will be shown subsequently, 
the observations of best comparison with data will change 
with Reynolds number. 

High Reynolds number Computations - Figures 17 through 
21 are high Reynolds number calculations showing the effect 
of turbulence model and viscous model on pressure coeffi- 
cient and turbulent kinetic energy distributions on the body. 
Figure 17 is a comparison of the four turbulence models at 
Nfa = 128 million using k-thin layer viscosity, min-mod lim- 
iter and tripl for turbulent tripping. The three ASM cluster 
around the experimental data matching the pressure recov- 
ery in the separated flow region considerably better than at 
low Reynolds numbers. The Standard K — e model predicts 
the strongest shock and highest pressure recovery. 

Figure 18 is the plot of peak turbulent kinetic energy 
similar to figure 13 for the four turbulence models. Signif- 
icantly. all four models developed turbulent flow immedi- 
ately downstream of the turbulent trip as seen by the four 
curves departing from the trip spike in K/a 2 at levels around 
0.004. Each turbulence model remained at slightly different 
levels, but had similar trends until the region of flow involv- 
ing the shock-separation downstream of x/d w = 0.25. The 
trend of the peak turbulent kinetic energy was similar to the 
7 million Reynolds number trend in figure 13. Though the 
three ASM have very similar static pressure coefficient dis- 
tributions, figure 17. the peak K/a 2 trends are completely 
different. Also, the Cp distributions between SZL and the 
Standard K — e model are very different, but the peak K/a 2 
have similar trends and levels. Therefore at this point, a cor- 
relation between the trend of K/a 2 and C p can not be made. 

Figure 19 is the effect of viscous model using the Girimaji 
ASM at 128 million Reynolds number. In this instance, 
the k-thin layer calculation (321) provides the best compar- 
ison with the experimental pressure coefficient distribution. 
The j-k viscous models behaved similarly in that the shock 
strength increased and the recovery pressure was higher than 
the k-thin layer calculation. Figure 20 is the peak K/ar for 
the three viscous models shown in figure 19. The three vis- 
cous models have similar trends in peak turbulent kinetic 
energy until the region of shock-separated flow downstream 
of x/d m = 0.25. Both j-k viscous models generate higher 
peak turbulence than the k-thin layer model. The plots in 
figures 21(a) to (c) are contours of turbulent kinetic energy 
predicted by the three viscous models previously discussed. 
The k-thin layer viscous model, figure 21(a). has an abrupt 
discontinuity in the flow-field around the boattail trailing 
edge. x/d m = .8. while the both j-k viscous models predict 
very smooth and continuous contours from the region of the 
shock. x/d m = .25. to downstream. 

Reynolds number Trends —Figures 22 through 27 are 
trends of integrated boattail pressure drag, skin friction, and 
predicted point of flow separation with Reynolds number. 
The integrated pressure drag variation with Reynolds num- 
ber comparing CFD with experiment is shown in figure 22. 
Despite the changes in the shock strength and pressures oil 
the nozzle boattail with Reynolds number: the variation in 
pressure drag was small. Overall, the predicted level of pres- 
sure drag was slightly below that of the experimental data, 
though at the low and high Reynolds numbers the CFD was 
almost within the scatter of the experimental data. As a 
point of reference. 3 additional data points are plotted to 
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include data obtained for the short cryogenic models tested 
in the 16 Foot Transonic Tunnel at Langley, and the origi- 
nal 48 inch model also test in the 16 Foot Tunnel. 

Figure 23 shows the predicted change in static pressure 
coefficient distribution with Reynolds number. The largest 
change seems to occur from the very low Reynolds number 
to the mid-range, with the code predicting a large increase 
in the peak velocity, a downstream shift of the peak and a 
slight elevation of the static pressure of the flow in the region 
of separation. Considerably less change was predicted be- 
tween the mid-range Reynolds number to the high Reynolds 
number of 128 million. 

Figure 24 is a bar chart of the integrated pressure drag 
oil the boattail at 7 million Reynolds number comparing 
the different viscous models and trip location predicted drag 
with the Girimaji ASM with experimental data. The higher 
recovery pressure that occurred through the j-k viscosity 
calculations reduced the integrated pressure drag from 37 
to roughly 28 nozzle drag counts. The scatter in the CFD 
results is about the same as the experimental results with 
the exception of the 48 inch model data tested in 16-Foot. 

High Reynolds number comparisons are shown in 
Figure 25 with the addition of GS, SZL and Standard K — e. 
The scatter in the CFD is similar to the low Reynolds num- 
ber comparison with the Standard K — e predicting the low- 
est drag due to the considerably higher pressure recovery at 
the boattail. Girimaji and SZL, k-thin layer, are the closest 
to the experimental data, though on the average are low. 

Variation of predicted skin friction coefficients for Girimaji 
ASM with Reynolds number is plotted against flat plate 
wetted area estimations in figure 26. In general, the CFD 
predicts skin friction coefficients are 3.5 nozzle drag counts 
low at 7 million Reynolds number and about 1.5 nozzle drag 
counts low at 128 million Reynolds number. Considering 
the flow effects not accounted for by the flat plate wet- 
ted area calculations, (e.g., non-constant Mach number, ad- 
verse/favorable pressure gradients, aft- projected areas and 
separated flow) this comparison is fairly good. 

Lastly, figure 27 is an analysis of the predicted point of 
flow separation with Reynolds number comparing with some 
flow visualization data obtained in the 16 Foot Transonic 
Tunnel on the 48 inch model in 1974 and the parametric 
theory of Reshotko-Tucker. ref. 18. The separation observed 
in ref. 18 was somewhat three dimensional with the esti- 
mated extent thereof shown by the spread in open triangles 
in the figure. No separation data are available for this model 
at any of the other Reynolds numbers. The SZL ASM pre- 
dicted a flow separation point that more closely matched the 
wind tunnel measurement and Reshotko-Tucker predictions 
with increasing Reynolds number. Both Girimaji and GS 
predicted flow separation points further downstream. The 
j-k viscosity predictions of Girimaji predicted the least sep- 
arated flows, with the j-k coupled viscosity calculation pre- 
dicting practically no separated flow at 7 million Reynolds 
number. 

Remarks 

1. The high Reynolds number boundary layer calculation of 

skin friction and shape factor for the subsonic flat plate 

was consistent with theoretically predicted behavior. 


2. The linear turbulence simulation predicted a shock fur- 
ther downstream and a recovery pressure higher than the 
non-linear turbulence simulations at the low and high 
Reynolds numbers. 

3. The best performance combination of turbulence mod- 
els and viscous models appears to change from low 
Reynolds number to very high Reynolds number. The 
ASM with j-k viscous modeling appeared to provide the 
best low Reynolds number comparison, while ASM with 
only k-thin layer viscosity most closely matched the high 
Reynolds number static pressure coefficient distribution. 
Further investigation is required to resolve this issue. 

4. The afterbody pressure drag variation observed in the 
experimental data and the computations with Reynolds 
number was small. The change with Reynolds number 
of the pressure coefficient distribution observed in the 
experimental data is qualitatively predicted by the CFD. 
This "no-effect effect" had been discussed in the previous 
high Reynolds numbers investigations. 

5. Most of the solutions using the non-linear models pre- 
dicted a separation point downstream of experimental 
flow visualization and parametric theory except the model 
by Shill, Zhu and Lumley. 
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Appendix 

The functions and variables used in the Girimaji algebraic 
Reynolds stress model are listed: 


1:I| =Cj + 2 


L,-*-*-. 

1 2 3 6 2 2 


m 


_ (KY q c . _(K 
V e / V e 


Wmi}W m n 


2L °i 


2 

(mi }) 2 


(L3) 2 + 2 t ?2 ( L4) 2 


•=(«- Yp-s( 2 ' , ’- 9M + 27r ) 

D = £ + £«*(«). 

4 27’ v/-aV27 


The coefficients G 2 and G 3 are: 
-L 4 G 1 


g 2 = 


Ll-r, x L\G{ Lq — t) X L\G\ 

additionally 

C x = 3.4: Cj = 1 . 8 ; C 2 = 0.36; C 3 = 1.25, C 4 = 0.4. 


2 L 3 Gi 
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M=0.097 

p w = 1 ,7229MPa 



|.e. *(m) 

Figure 1 .-Block and grid arrangement for subsonic flat plate. 



x(m) 

Figure 2.-Reynolds number with distance along flat plate. 



y + y* 

(a) Law-of-the-wall profile at R e =1420. (b) Turbulent shear stress at R e =1420. 

Figure 3 - Comparison of boundary layer characteristics for different turbulence models. 


10 


American Institute of Aeronautics and Astronautics 



Girimaji ASM 
Gatski-Speziale ASM 
Shih.Zhu & Lumley 
Laminar Sublayer 
Buffer Region 
Turbulent Region 
Turbulent Region-2 

jid 


(a) Law-of-the-wall profile at R 0 =1 00,000. 



(b) Turbulent shear stress at R e =1 00,000. 


Figure 4.- Comparison of boundary layer characteristics for different turbulence models. 
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Figure 5. -Local skin friction for subsonic flat plate. 



Figure 6. -Average skin friction for subsonic flat plate, (same symbol table as figure 5). 
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Figure 7.- Turbulent viscosity development 
with R 61 
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Figure 8.- Boundary layer shape factor development 
with R x 



Figure 9.- Photograph of 8 inch model in 0.3m Transonic Cryogenic Tunnel. 
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Figure 10.- Representative grid sensitivity at 
7 million Reynolds number. 
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Figure 1 1 Representative grid sensitivity at 
128 million Reynolds number. 
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N Re =7mil.,Std.K-E,321,trip1 

N Re =7mil.,Girimaji,321,trip1 

N R .=7mil.,GS,321,trip1 

N He =7mil. 1 SZL,321,trip1 

O L=1 6,M=0.90,7.028M,4>=0. 

D L=16,M=0.90,7.028M,4»=120. 

A L=1 6,M=0.90,7.028M,<J>=240. 

• L=16,M=0.90,7.107M,4>=0. 

■ L=1 6,M=0.90,7.1 07M,<t>=1 20. 

A L=16,M=0.90,7.107M,<t)=240. 


Figure 12.- Comparison of turbulence models with experimental data, M=0.9,N Re =7million. 
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Figure 16.- Effect of viscous model on nozzle 
C p distribution, N Re = 7 million. 



x/d m 


Figure 17.- Comparison of turbulence models 
with data, N Re =128 million. 



Figure 18.- Peak turbulent kinetic energy in boundary layer, M=0.9, N Re =128 million. 
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Figure 1 9.- Effect of viscous model on nozzle C p distribution, N R6 =1 28 million. 
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Figure 20.- Effect of viscous model on peak turbulent kinetic engery, M=0.9, N Re =128 million. 
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(c) j-k uncoupled viscous model. 

Figure 21 Effect of viscous model on turbulent kinetic energy contours. 
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Figure 22.- Comparison of predicted integrated pressure boattail drag with experiment, M=0.9. 
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Figure 23.- Predicted variation of afterbody C p distribution with Reynolds number, M = 0.9. 
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Figure 24.- Comparison of integrated boattail pressure drag at around 7 million Reynolds number, 
(PAB3D results using Girimaji ASM), M =0.9. 
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Figure 25.- Comparison of integrated boattail pressure drag near 128 million Reynolds number, 
M =0.9. 
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Rgure 26.- Comparison of predicted skin friction coefficient with wetted-area, flat-plate estimation, M = 0.9. 
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Figure 27.- Variation of predicted point of flow separation compared with flow visualized data and 
parametric theory. 
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